

#include "glodef.h"
#include "errlog.h"
#include "quadr.h"

#include <stdlib.h>
#include <string.h>

#define QD2     1
#define QD4     2
#define QD5     3

const int QDEG2=QD2;
const int QDEG4=QD4;
const int QDEG5=QD5;

/* The following formula is degree 2 */
/*-----------------------------------------------------------------------*/
/* the quadrature nodes, weights and numbers of them */
static const int QN_QD2=5;
static const double QPts_QD2[5][3]={
    {0.0, 1.0/3.0, 1.0/3.0},
    {1.0/3.0, 0.0, 1.0/3.0},
    {1.0/3.0, 1.0/3.0, 0.0},
    {1.0/3.0, 1.0/3.0, 1.0/3.0},
    {0.25, 0.25, 0.25}
};
static const double QWts_QD2[5]={3.0/40.0, 3.0/40.0, 3.0/40.0, 3.0/40.0, -2.0/15.0};
/*-----------------------------------------------------------------------*/


/* The following formula is degree 4 */
/*-----------------------------------------------------------------------*/
/* the quadrature nodes, weights and numbers of them */
static const int QN_QD4=11;
#define A2      0.399403576166799219
#define B2      0.5-A2
#define W2     -0.131555555555555555e-1
#define W2A     0.762222222222222222e-2
#define W2B     0.248888888888888888e-2
static const double QPts_QD4[11][3]={
    {0.25, 0.25, 0.25},
    {1/14, 1/14, 1/14},
    {11/14, 1/14, 1/14},
    {1/14, 11/14, 1/14},
    {1/14, 1/14, 11/14},
    {A2, A2, B2},
    {A2, B2, A2},
    {B2, A2, A2},
    {B2, B2, A2},
    {B2, A2, B2},
    {A2, B2, B2}
};
static const double QWts_QD4[11]={W2, W2A, W2A, W2A, W2A, W2B, W2B, W2B, W2B, W2B, W2B};
/*-----------------------------------------------------------------------*/

/* The following formula is degree 5 */
/*-----------------------------------------------------------------------*/
/* the quadrature nodes, weights and numbers of them */
static const int QN_QD5=15;
#define A3      0.433449846426335728
#define B3      0.665501535736642813e-1
#define W3A     0.602678571428571597e-2
#define W3B     0.302836780970891856e-1
#define W3C     0.116452490860289742e-1
#define W3D     0.109491415613864534e-1
static const double QPts_QD5[15][3]={
    {0.0, 1.0/3.0, 1.0/3.0},
    {1.0/3.0, 0.0, 1.0/3.0},
    {1.0/3.0, 1.0/3.0, 0.0},
    {1.0/3.0, 1.0/3.0, 1.0/3.0},
    {0.25, 0.25, 0.25},
    {8.0/11.0, 1.0/11.0, 1.0/11.0},
    {1.0/11.0, 8.0/11.0, 1.0/11.0},
    {1.0/11.0, 1.0/11.0, 8.0/11.0},
    {1.0/11.0, 1.0/11.0, 1.0/11.0},
    {A3, A3, B3},
    {A3, B3, A3},
    {B3, A3, A3},
    {B3, B3, A3},
    {B3, A3, B3},
    {A3, B3, B3}
};
static const double QWts_QD5[15]={W3A, W3A, W3A, W3A, W3B, W3C, W3C, W3C, W3C, W3D, W3D, W3D, W3D, W3D, W3D};

/*-----------------------------------------------------------------------*/


void integrate_ref_el(int Q, int N, __func_t F, double *V, void *data)
{
    int i, j;
    double *R;
    int QN;
    const double (*QPts)[3], *QWts;
    switch(Q)
    {
        case QD2: QN=QN_QD2;
            QPts=QPts_QD2;
            QWts=QWts_QD2;
            break;
        case QD4: QN=QN_QD4;
            QPts=QPts_QD4;
            QWts=QWts_QD4;
            break;
        case QD5: QN=QN_QD5;
            QPts=QPts_QD5;
            QWts=QWts_QD5;
            break;
        default: error_msg("Quadrature formula %d not implemented.\n", Q);
    }

    /* clean up memory */
    memset(V, 0, N*sizeof(double));

    /* compute the quadrature formula */
    for(i=0; i<QN; i++)
    {
        F(QPts[i], &R, data);
        for(j=0; j<N; j++)
            V[j]+=QWts[i]*R[j];
    }

}

